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Macroscale experiments show that a train of two immiscible liquid drops, a bislug, can spon- 
taneously move in a capillary tube because of surface tension asymmetries. We use molecular 
dynamics simulation of Lennard- Jones fluids to demonstrate this phenomenon for NVT ensembles 
in sub-micron tubes. We deliberately tune the strength of intermolecular forces and control the 
velocity of bislug in different wetting and viscosity conditions. We compute the velocity profile of 
particles across the tube, and explain the origin of deviations from the classical parabolae. We show 
that the self-generated molecular flow resembles the Poiseuille law when the ratio of the tube radius 
to its length is less than a critical value. 



I. INTRODUCTION 

Dynamics of liquid flow in nanoscopic systems has 
received great attention in recent years. Design and 
manufacturing of nanofluidic devices, fluid transport 
phenomenon^, the recent development of "lab-on-a- 
chip" technology 3 ^ and fluid imbibition in the nanopores 
of biomembranes, are amongst notable applications. Our 
knowledge, however, has been mostly confined to the dy- 
namics of capillary flows at macroscales^ where the re- 
lationships between the velocity and contact angle with 
the surface tension, dynamic viscosity and wetting char- 
acteristics are known 7 .^. 

Generating and controlling the driving force of liq- 
uids in extremely thin tubes is a major technologi- 
cal challenge, which suggests the application of self- 
propelling systems. The spontaneous motion of liquid 
drops was initially reported by Marangoni^ and it was 
shown to occur in variably wet and temperature con- 
trolled tube o 10 ^ 1 , chemically reacting fluid s 12 ' 13 ? 14 and in 
surface-tension-driven motions^. In all these cases the 
driving force comes from physical or environmental asym- 
metries, which may cause a net motion. In the present 
study, we are interested in the surface-tension-driven mo- 
tion of droplets because of its importance in microflu- 
idic automation 8 -^. In particular, we intend to scale 
down (by several orders of magnitude) a self-propelling 
liquid system, called a bislug, which was introduced by 
Bico & Quer o 15 ' 17 . Although the motion of a bislug is 
fairly explained by the dynamics of Newtonian fluids in 
macroscales, there is neither theoretical nor experimental 
evidence for the applicability /validity of the same phys- 
ical rules for molecular ensembles at sub-micron scales. 

A liquid bislug is composed of two liquid droplets, liq- 
uid I and liquid II, juxtaposed inside a capillary tube and 
bounded by a third fluid (fluid III) on both sides (Fig. 
[T]a) . The system therefore involves three interfaces each 
called a meniscus. The difference between adhesive inter- 
actions of two neighboring fluids directly determines the 
meniscus wetting power, and so its contact angle. In the 
partial wetting state, the contact angle is finite, and it re- 
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FIG. 1: (a) Configuration of a liquid bislug in a nanotube. 
(b) Initial state of our LJ fluids. Fluid III initially fills the 
tube and liquids I and II occupy the reservoir. 



duces to zero in the limit of complete wetting condition. 
Moreover, the cohesive interactions inside each of two 
adjacent fluids as well as adhesive interactions between 
them, determine the meniscus surface tension. The total 
driving force exerted on a bislug is 

F cx 7 d — 7i cos 6»i + 712 cos 6 12 - 72 cos 2 , (1) 

where 7^ is the effective surface tension, and the param- 
eters (71,712,72) and (0i,0i2,02) are, respectively, the 
surface tensions and contact angles (Fig. [TJa). 

In macroscale experiments, when all menisci partially 
wet the tube wall, the contact angles are changed so that 
the forces at three menisci compensate each other and 
keep the liquid bislug in static equilibrium&i 7 -. As soon 
as at least one of the menisci completely wets the tube, its 
contact angle drops to zero and the existing force balance 
breaks down. If liquid II is more wetting than liquid 
I (see Fig. [T]a), it will be deposited on the tube wall 
while liquid I moves with a constant velocity. In this 
article, we carry out a molecular dynamics simulation of 
a nanoscopic system of immiscible liquids, and reveal its 
self-propelling dynamics for a wide range of material and 
intcrfacial properties. 
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II. MODEL DESCRIPTION 

We use the setup of Fig. [IJ) to incrementally turn on 
the force contribution of each meniscus. Using fluid III, 
we simulate the presence of air or any other fluid and 
equalize its pressure at both sides of the bislug using 
periodic boundary conditions. Fluid III also suppresses 
possible diffusion effects. We consider a cylindrical tube 
of the length L and radius R composed of fixed solid par- 
ticles that do not interact with each other. The wall par- 
ticles are distributed on a cylindrical arrangement with 
a constant distance of a, which is our length scale. This 
is how wc build a rigid impenetrable wall. The tube is 
attached to a (L x — L) x L y x L z (L y = L z = 2.5R) boxy 
reservoir on the left side. The height (and also width) 
of the reservoir is larger than the tube diameter, and the 
joint between the tube and the reservoir is rigid. Wc 
apply periodic boundary conditions at the lateral walls 
of the reservoir (perpendicular to the view shown in Fig. 
QJ>) as well as in the longitudinal direction: all liquid par- 
ticles may freely enter the reservoir if they reach to the 
rightmost boundary of the tube. 

There are four different types of particles correspond- 
ing to liquid I, liquid II, fluid III, and the solid wall, which 
are numbered in our subsequent formulation as 1, 2, 3, 
and 4, respectively. Fluid particles are located initially 
at a cubic lattice with the lattice constant of er, and all 
particles interact through the Lennard-Jones (LJ) poten- 
tial 



U ij (r)=Ae i 



(2) 



Here and define the strength and effective range 
of Uij between the particles of types i and j. Theoreti- 
cally, the radius of influence of the LJ potential extends 
to infinity and all binary interactions will involve non- 
zero finite forces. Nevertheless, the computational cost 
will be drastically decreased if one cuts off the attractive 
LJ force for radii above some r c . For all interactions, we 
choose r c — 2.5er, which has already been used in the 
literatur o 18 i 19 ' 20 and it guarantees a drop of the attrac- 
tive force below 2% of its maximum value. The choice 
of r c — 2.5a has also resulted in a reasonable match be- 
tween molecular dynamics simulations and experimen- 
tal data^i. We control the mixing of different fluids by 
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A; 



1 if i = j, 
0.1 if Mi. 



(3) 



for i,j = 1,2,3. 

The surface tension at each interface depends on the 
internal interactions of each fluid as well as interactions 
between two fluids, controlled by e.y = 1,2,3). The 
wetting power of each liquid depends on adhesive inter- 
actions between the liquid and tube particles, which is 
solely controlled by ea in the LJ potential. The dif- 
ference between the values of en/en for adjacent fluids, 



determines the contact angle of their interface (menis- 
cus) such that the concave part faces towards the fluid 
with lower en/ en. The resultant driving force F is con- 
trolled by adjusting the values of and so the contact 
angle a 19 ' 20 . Moreover, the viscosity of the ith fluid, rji, is 
determined by en. 

We implement non-equilibrium molecular dynamics 
simulation of an NVT ensemble^ 3 , using the velocity Ver- 
lct algorithm with the integration time step St = 0.005io 
where 



to = 



1 



48fc B T V38X 



(4) 



The chosen time step corresponds to 10~ 14 seconds for a 
wide range of fluids, and its efficiency has already been 
confirmed in the study of liquid argo n 23 ' 24 . We use the 
dimensionless values of a = 1, the particle mass m = 1, 
and set k^,T = 0.8. To preserve the ensemble tempera- 
ture at a constant value, a Nose-Hoover thermostat 2 ^ is 
applied at each step. We scale all lengths and velocities 
by a and cr/to, respectively. 



III. SIMULATION RESULTS 

We do our simulations when: (A) all fluids have the 
same viscosity 77 and surface tension 7 (B) viscosities of 
fluids differ from each other. We then generate the nor- 
malized displacement x/L x of the front meniscus (FM) 
versus the normalized time t/ta for different values of 
a = R/L. The position of the front meniscus x is mea- 
sured from the leftmost boundary of the reservoir, and 
it is defined as the location where the density of liquid 
I reaches to 90% of its maximum value. Although the 
experiments of type A are unrealistic, they help us to 
independently investigate the influence of contact angles 
on the bislug dynamics. We study three possible circum- 
stances: (i) only the rear meniscus (between liquid II 
and fluid III) completely wets the tube (ii) the rear and 
middle menisci completely wet the tube (iii) all three 
menisci completely wet the tube. We use two tubes of 
radii R = 10.34 and 5.17. 



A. Equal Viscosities 

We set eii=e22=£33=l to create equal viscosities, and 
use £34 = 0.5 and ei2=£23=£i3=l to get the same bound- 
ary interactions between different fluids. The internal 
interactions of each fluid is governed by /3u=l (i=l,2,3) 
and we set /3i2=/323=/?i3=0.1— to prevent liquids from 
mixing. We start our computations with a case that all 
menisci partially wet the tube. The contact angles in a 
partial wetting condition lie in the range < 6\, 6*i 2 , 9 2 < 
90° should one choose the values of 0.5 < ei4,e 2 4 < 1.1. 
In such a circumstance, the capillary force pulls both liq- 
uids into the tube subsequently, while fluid III leaves the 
tube and enters the reservoir. 



3 



0.8 


- 8,4=0.5 
8,4=0.6 

- ^ 8 24 =0.7 






K 

0.6 

\ 

*\ 


^_ 8,4=0.9 

□ e 24 = 1.3 /2 

8 24 = 1.4 // 

- — £24=1.5 J J 






0.4 












8 14 =0.7 




0.2 





i 






2000 4000 6000 8000 

t/to 

FIG. 2: The normalized displacement curve of the front 
meniscus (x/L x ) versus time in partial and complete wet- 
ting conditions with a = 0.107(7? = 10.34, L = 96.14) and 
T]i = 772 = 773. 

As each meniscus enters the tube, the driving force is 
altered and the slope of the displacement curve of the 
FM varies accordingly. One can therefore identify three 
phases on displacement curves (Fig. The first phase 
begins when liquid I enters the tube and it ends when 
liquid II flows into the tube. The driving force of the 
first phase is proportional to jd — 7COS#i. The sec- 
ond phase starts by the inflow of liquid II and contin- 
ues until the last molecules of liquid II leave the reser- 
voir. This stage involves two active menisci implying 
Id = 7(cos#i + cos#i 2 ). In the third (and last) phase, 
the whole bislug is inside the tube and the driving force 
is computed from Eq. JT]). 

Fig. [2] displays the displacement history of the FM for 
several values of 0.5 < £24 < 0.9 and for £14 = 0.7. By 
keeping £14 constant, we have fixed the driving force of 
the FM. For the three cases £24=0.5, 0.6, and 0.7 that 
satisfy £24 < £14 = 0.7, the concave side of the mid- 
dle meniscus (between liquids I and II) faces liquid II 
(#12 > 90°) and it resists the driving force of FM. Con- 
sequently, in the second phase the bislug moves slower 
than the first one. For £24 = 0.5 we have £24= £34, 
whose natural implication is £14 — £24= £14 — £34, which 
in turn leads to the condition 6*12 + 6\ = 180°. This 
says that the surface tension of the front meniscus can- 
cels that of the middle one, and the bislug must stop 
as soon as liquid I leaves the reservoir (zero velocity in 
the second part of displacement curve in Fig. HJa). By 
increasing £24, the velocity of FM increases as long as 
liquid II has not completed its entrance into the tube. 
However, after the complete drainage of liquid II from 
the reservoir, the displacement curve of the FM becomes 
flat indicating a static equilibrium in partial wetting con- 
ditions: 7(cos6>i + cos6>i2 — cos6> 2 ) = 0. 

For £24 > 1.2 and with the same £14 = 0.7, the rear 
meniscus completely wets the tube, and results in 82 = 
and F oc 7^ = 7(cos6*i + cos #12 — 1). This is the first 
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FIG. 3: The normalized displacement curve of the front 
meniscus (x/L x ) versus time in partial and complete wet- 
ting conditions with a = 0.107 (-R=10.34, L=96.14) and 
771 = 772 = 773. 




FIG. 4: A moving bislug corresponding to the case of ei4=0.7 
and £24 = 1-5 in Fig. [2] 



complete wetting situation, which holds over the range 
1.2 < £24 < 1.5. Since cos#i + cos6*i 2 > 1, the bislug 
undergoes a sustained motion after entering the tube un- 
til liquid II is totally consumed, or the FM reaches the 
end of the tube (Fig. [2]). By further increasing £24, we 
obtain a critical value of £24 ~ 1.5 beyond which the 
middle meniscus completely wets the tube. This is the 
second complete wetting state (82 = 6*12 = 0) and the 
force components corresponding to 7 cos 82 and 7 cos 812 
cancel each other identically. The bislug is thus driven 
by the FM: F oc 7^ = 7Cos6>i. Fig. 2] shows a snapshot 
of the moving bislug in its third phase. It is evident that 
liquid II leaves a thin stationary layer on the tube wall. 

To obtain the state of complete wetting condition for 
all three menisci, we fix £14 = 1.2 and carry out new 
simulations (see Fig. [3]). For the very special case of 
£14 = £24 = 1-2 that implies 0\ = 8 2 = and 6>i 2 = 90°, 
and according to our assumption 7=71=712=72; liquids 
I and II will behave as a single droplet, which is pulled 
into the tube by the maximum force F\ cx 7 exerted by 
the FM, but it stops once the rear meniscus enters the 
tube and contributes Fn oc —7 to the total driving force. 
We then increase £24 so that the force components 7 cos 8\ 
and 7 cos 82 cancel each other and the driving force comes 
from the middle meniscus whose concave side faces liq- 
uid I for £24 > £14. The middle meniscus, however, will 
be in partial wetting as long as cos #12 < 1. Therefore, 
the velocity of the FM (and so the bislug) is boosted by 
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increasing £24 because cos #12 oc £24 — £14- The four last 
graphs of Fig. [3] show that the velocity of the FM is sat- 
urated for £24 > 2.8, which corresponds to the complete 
wetting condition of all menisci. 

The experiments with (£14, £24) = 0.7 in Fig. [2] and 
(£14,624) = 1.2 in Fig. [3] show that the bislug velocity 
in the second phase is slightly less than the first phase 
while the equality of £14 and £24 should have resulted 
in 6*12 = 90° and identical velocities during the first 
and second phases of motion. Our explanation of this 
phenomenon is as follows: while the middle meniscus is 
driven, it is subject to back-and-forth oscillations that al- 
ternately change the contact angle and keep its average 
magnitude larger than the expected value (this is because 
the particles near the contact point of a meniscus usu- 
ally lag the field particles, which feel less frictional force). 
For the special case of £14 = £24, the average value of #12 
marginally exceeds 90° and causes a tiny resistive force. 
Therefore, the velocity of the bislug is slightly reduced 
during the second phase. So far we have considered dif- 
ferent cases of bislug motion when three fluids have the 
same viscosity 77. Below we vary the viscosities of fluids 
and investigate realistic cases. 

B. Different Viscosities 

We use Kirkwood-Buff 's expression 

7 = j dxy>xx(x) ~ ^ [Pyy(x) + Pzz{x)]j , (5) 

and compute the surface tension 7 from the elements of 
pressure tensor using the equilibrium molecular dynam- 
ics simulation of two adjacent, thermostated liquids in a 
cube2£. In similar conditions, we apply the Green-Kubo 
integral to obtain the viscosity of our LJ fluid a 23 i 27 . The 
element p X x{x) of the pressure tensor is normal to the 
interface. For the system of fluids with equal viscosi- 
ties and surface tensions, we find 7 = 2.6±0.15 and 
77 = 3.25 ± 0.25. It is remarked that the application 
of other methods such as shear flow, gives compatible 
results^. 

We now vary the viscosities of fluids around 77 w 3.25. 
Increasing the value of en enables us to increase the vis- 
cosity of the ith fluid, rji. By taking £34 = 0.1 and 
e«=(0.5, 0.8, 1, 1.1, 1.2), we find 

rji = (1.6±0.1,2.6±0.2,3.25±0.25,3.6±0.32,4.±0.33), (6) 

and apply our NVT code to compute x(t)/L x when 
771 > 3.25 > r)2 > ^3 = 1.6. Some of our results arc 
displayed in Fig. [5l As one could anticipate, the first 
two phases of displacement curves are no longer linear 
versus time, and the profile of the first phase resembles 
a parabolic pattern x/L x oc y/t/to, which is the well- 
known Washburn equation^. It is evident that even for 
the case of different viscosities, the third phase of each 
displacement curve has still a constant velocity in accor- 
dance with macroscale experiments^. 
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FIG. 5: The normalized displacement curve of the front 
meniscus (x/L x ) versus time for a = 0.107 (i?=10.34, 
L=96.14), (ei4, £24) = (0.7,2.3) and different viscosities 771 > 
T]2 > 773- Here, the FM is in partial wetting condition and the 
other two menisci completely wet the tube wall. We demon- 
strate the special case of 771=772 > 773 to show the effect of 
771 and 772 on the bislug velocity independently. Horizontal 
dashed line approximately indicates the beginning of the third 
phase. 

By comparing the results of Figures [2] and [5] we con- 
clude that reducing 773 (with respect to 771 and 772) boosts 
the bislug velocity considerably. It is remarked that we 
needed to decrease £34 while we were doing so for £33 
just to preserve the wetting condition of the FM. In 
Fig. El comparing the graphs of (fin, £22) = (1,0.8) and 
(eilj £22) = (1, 1) suggests that decreasing 772 (while 771 is 
kept constant) induces an even faster motion because the 
deposition of liquid II on the tube wall speeds up. The 
graphs corresponding to £22 = 0.8 show that the bislug 
velocity is slowed down by increasing en, which means 
higher values of 77! and therefore a larger viscous force. 
This effect is clearly evident from the longer first phase 
(similar to \ftjta) of the bislug displacement curve. 

IV. PHYSICS OF THE FLOW 

In all cases, the velocity of the FM, V(t) = dx/dt, is 
constant over the last phase of displacement curves, and 
it is proportional to jd- The motion with a constant V is 
due to the balance between the driving force F = 2TrR-fd 
and the viscous friction force Fdrag from the tube wall, 
i.e., F — Fdrag = mdV/dt = where m is the total mass 
inside the tube. When all fluids have the same viscosity 
77 and the flow is Poiseuille, one obtains 

V = V p = ^, F diag = 8nr]LV p . (7) 

In this equation, 77 and L must be replaced, respectively, 
with 77! = 772 and the bislug length when fluid III has a 
negligible contribution to the dynamics. Given a = R/ L, 
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the magnitude of V p can therefore be calculated. To see 
in which conditions the bislug velocity attains the value of 
V p , we carry out some new calculations with the two radii 
R = 5.17 and 10.34, and by taking several tube lengths 
of 52.51 < L < 96.14, which yield 0.107 < a < 0.197. 
We use fluids of the same viscosity, set £14 = 1.2, and 
let all menisci completely wet the tube wall. Some of 
our results have been illustrated in Fig. [5J According 
to our simulations, the quantity V/L x tends to V p /L x 
by decreasing a and the relative error \V — V p \/V p re- 
duces from 33% (for a = 0.197) to 3% (for a = 0.107). 
In the complete wetting condition of all menisci and for 
£24 = 3.5, we find the critical value a cr ~ 0.1 below 
which the velocity of the FM becomes identical to V p . A 
couple of our experiments correspond to R — 5.17 and 
L = 96.14 that give a = 0.053. With such a small ratio 
and for £14 = £24 = 1.2, liquids I and II remain sepa- 
rated but they behave as a single drop. For the stronger 
wall interactions of liquid II with £24 = 3.5, the relative 
velocity error further decreases and we obtain V ~ V p . 
This is evident from Fig. [5] as the first and third phases 
of motion have the same slope. 

We are using a Nose-Hoover thermostat, which is 
mainly applied to equilibrium phenomena. The self- 
propelling motion of liquid trains, however, is a non- 
equilibrium case and a thermostat may destroy the 
Galilean invariance of dynamical equations. Through 
computing the Washburn prcfactor, we show that such 
an effect is quite small in most of our simulations that 
correspond to narrow tubes of a < a cr . When fluids I 
and III are immiscible and have the same viscosity, the 
time dependence of the displacement function x/L x is lin- 
ear during the entrance of liquid I into the tube, and the 
Washburn equation V ~ < -1 / 2 is transformed to the ordi- 
nary law of capillary dynamics: V = 7a cos 6/ (47/). This 
becomes V = 0.2a for (7, ?y)=(2.6, 3.25) in our setup and 
in the complete wetting condition of the FM (9 = Q\ =0). 
For £14 = 1.2 and a = 0.107, the first phase of the 
displacement curves in Fig. [3] corresponds to the com- 
plete wetting of the FM and the slope of this phase gives 

V w 0.0225. On the other hand, the capillary equation 
yields V = 0.2 x 0.107 = 0.0214. For the case with 
(R,a) = (5.17,0.053) and e lA = £ 24 = 1.2 in Fig. © the 
capillary equation leads to V = 0.2x0.053 = 0.0106 while 
the slope of the first phase of the displacement curve is 

V fw 0.0108. The relative errors between theoretical pre- 
dictions and our simulations are therefore 5% and 1.8% 
for a = 0.107 and a = 0.053, respectively. This result 
shows that the effect of the Nose-Hoover thermostat on 
the Galilean invariance of governing equations is insignif- 
icant, specially when a < a cr . 

Let v((, r, t) denote the longitudinal velocity of 
molecules inside the tube where r and £ are the distances 
of molecules from the tube centerline and the reservoir 
outlet, respectively, and t is the time. Small circles and 
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FIG. 6: The normalized displacement curve of the front 
meniscus (x/Lx) versus t/to for 0.107 < a < 0.197 with 
R = 10.34, and for a = 0.053 with R = 5.17. All menisci 
are in complete wetting condition, and we have set 614 = 1.2 
and 7/1 = 772 > 773 • 

squares in Fig. [7] display the simulated profile of 

v(r) = lim - [ dty [ * «(C, r, t)dC, (8) 

together with the best-fitted parabolic curves to the 
data. These data have been obtained for the case with 
(£i4,£ 24 ) = (0.7,1.5) in Fig. S and (e u , e 22 )=(l, 0.8) in 
Fig. [5l The time average in equation (jHJ) is computed 
only over the third phase of bislug motion, and averag- 
ing along the tube length helps us to reduce the noise of 
dynamic capillary waves and work with a better statisti- 
cal data. It is seen that for a = 0.107, the flow inside the 
tube is almost Newtonian and laminar. There is a bet- 
ter match between the simulated data and the best-fitted 
parabolic curve near the wall because the boundary layer 
(due to shearing) initially develops near the wall and then 
penetrates into inner regions. In other words, molecules 
near the tube centerline are subject to radial motions 
more than those living near the wall. Our numerical ex- 
periments show that the functional form of v(r) does not 
depend on the contact angle 812 but on a. By increasing 
a, the gradient 4-v(r) increases near the wall, enhances 
Fdrag and we get V < V p . The existence of a cr has also 
been confirmed by the macroscale experiments of Bico & 
Queroii. They report a cr rts 0.07 and ascribe the drop 
of velocity for a > a CT to the viscous dissipation at the 
menisci. Below, we argue that the physical mechanism 
of the process is different in molecular ensembles. 

The information of any developing phenomenon, like 
meniscus oscillation and boundary layer formation, is 
transmitted by individual molecules that interact by their 
neighbors and after several "kicks" lose their initial mem- 
ory. The oscillations of driving menisci remain unseen 
for downstream particles when the tube length is in- 
creased: for a long bislug, the number of longitudinal in- 
teractions between molecules is quite large and the com- 
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FIG. 7: The profile of v(r) versus the radial distance r from 
the tube centerline. Circles correspond to the case of equal 
viscosities with (£14,624) = (0.7, 1.5) in Fig. [2] Squares show 
the profile of the case with (en, £22) = (1,0.8) in Fig. [5] (for 
different viscosities). Solid and dashed lines are the best-fitted 
parabolic curves. 

ing information from oscillating menisci is gradually dis- 
solved. What actually matters in the radial direction 
is the development of the boundary layer whose propa- 
gation into central regions is again through interactions 
between molecules. When the tube radius is small, the 
number of radially disturbing kicks that molecules expe- 
rience is small and the boundary layer can fully extend 
to the tube centerline. In tubes of larger radii, however, 
molecules undergo too many radial interactions (kicks) 
and lose the "shearing information" that they had to 
carry to molecules near the tube centerline. As a re- 
sult, the boundary layer is not fully developed and some 
amount of the kinetic energy is carried by radial motions. 
An evidence for this process is the flattening of the pro- 
file of v(r) in central regions (see Fig. [7]). To have a 
Poiseuille flow in sub-micron scales, we thus need both 
a small radius (measured in units of a) and a large tube 
length. These would mean the existence of some a cr , 
which should correlate with the relaxation time of the 
bislug. 



V. CONCLUSIONS 

The bislug motion is a consequence of the transfer of 
free energy, from initially separated liquids, to the lon- 
gitudinal component of the kinetic energy of molecules. 
The system chooses a definite direction of motion because 
the different wetting capabilities of liquids I and II breaks 
the symmetry along the tube axis. Independent of the 
viscosities of fluids, we observed a constant bislug veloc- 
ity when both liquids I and II leave the reservoir. The 
Washburn law governs the history of x(t/to)/L x in the 
entrance phase of liquid I to the tube^. The bislug will 
not move (the first four cases in Fig. \2§ when the adhe- 
sion with the wall is not enough for making a permanent 
and stable binding between all molecules of a component 
(cither of liquids I and II) and the tube wall. 

We showed that the spontaneous motion of immiscible 
liquids inside the tubes of sub-micron radii is an efficient 
process in certain wetting conditions, and it leads to sig- 
nificant mass transport should one use appropriate liq- 
uids. An external pressure gradient has played a decisive 
role in previous experiments^ and simulations^ of water 
flow in carbon nanotubes. Our fundamental achievement 
was to capture the flow characteristics of a nanoscopic 
self-propelling system in the absence of an external pres- 
sure field. The size of the bislug, which indeed determines 
the relaxation time, and its relation with the form of the 
velocity profile across a small tube, helped us to suggest 
a new physical mechanism for explaining the deviations 
from the Poiseuille law. Our results have potential ap- 
plications in mass transport in fluidic devices and also in 
the separation of immiscible liquids in sub-micron scales. 



VI. ACKNOWLEDGMENTS 

We thank the referees for their constructive comments, 
which helped us to improve the presentation of our sim- 
ulations. This work was partially supported by the re- 
search vice-presidency at Sharif University of Technology 



* Electronic address: mjalali@sharif.edu 

1 M. Whitby and N. Quirke, Nature Nanotechnology 2, 87 
(2007). 

2 S.W.P. Turner, M. Cabodi, and H.G. Craighead, Phys. 
Rev. Lett. 88, 128103 (2002). 

3 H. Craighead, Nature 442, 387 (2006). 

4 D. Mijatovic, J.C.T. Eijkel, and A. van den Berg, Lab on 
a Chip 5, 492 (2005). 

5 E.W. Washburn, Phys. Rev. 17, 273 (1921). 

6 L.J. Yang, T.J. Yao, and Y. Tai, J. Micromech. Microeng. 
14, 220 (2004). 

7 PG. de Gennes, Rev. Mod. Phys. 57, 827 (1985). 

8 W.K. Chan, and C. Yang, J. Micromech. Microeng. 15, 



1722 (2005). 

9 C. Marangoni, Ann. Phys. Chem. 143, 337 (1871). 

M.M. Weislogel, AIChE J. 43, 645 (1997). 

1 V. Mazouchi and G.M. Homsy, Phys. Fluids 12, 542 
(2000). 

2 PG. de Gennes, Physica A 249, 196 (1998). 

3 C. Bain, G. Burnett-Hall, and R. Montgomerie, Nature 
372, 414 (1994). 

4 F. Domingues dos Santos and T. Ondarcuh, Phys. Rev. 
Lett. 75, 2972 (1995). 

5 J. Bico and D. Quere, Europhys. Lett. 51, 546 (2000). 

6 T.M. Squires and S.R. Quake, Rev. Mod. Phys. 77, 977 
(2005). 



7 



J. Bico and D. Quere, J. Fluid Mech. 467, 101 (2002). 
G. Matric, F. Center, D. Seveno, D. Coulon and J. De 
Coninck, Langmuir. 18, 7971 (2002). 

D.I. Dimitrov, A. Milchev, and K. Binder, Phys. Rev. Lett. 
99, 054501 (2007). 

D. I. Dimitrov, A. Milchev, and K. Binder, Langmuir 24, 
1232 (2008). 

M. Rao, and D. Levesque, J. Chem. Phys. 65, 3233 (1976). 

E. Diaz-Herrera, J. Alejandre, J. Ramirez- Santiago, and 

F. Forstmann, J. Chem. Phys. 110, 8084 (1999). 

M.P. Allen, and D.J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, Oxford, 1987). 
D. Frenkle and B. Smit, Understanding Molecular Simu- 
lation, From Algorithm to Applications (Academic Press, 



2002). 

W.G. Hoover, Phys. Rev. A 31, 1695 (1985). 

J.G. Kirkwood and F.P. Buff, J. Chem. Phys. 17, 338 

(1949). 

J.S. Rowlinson and B. Widom, Molecular Theory of Cap- 
illarity (Clarendon, Oxford, 1982). 

J. A. Backer, CP. Lowen, H.C.J. Hoefsloot, and P.D. 
Ledema, J. Chem. Phys. 122, 154503 (2005). 
M. Majumder, N. Chopral, R. Andrews, and B.J. Hinds, 
Nature 438, 44 (2005) 

J. A. Thomas and A.J.H. McGaughey, Phys. Rev. Lett. 
102, 184502 (2009). 



